Journal of the Royal Society Interface

With the development of microfluidics, it has now become possible to assess the susceptibility of bacteria to antibiotics at the single-cell level instead of relying on population measurements. Such studies are particularly relevant when the growth of bacterial population in the presence of antibiotics is heterogeneous. Here, we build a model to describe such a case, and apply it to experimental measurements on a small population of E. Coli exposed to ciprofloxacin, a drug which is well known for triggering a bistable response.

Time-resolved measurements are central to calibrating mechanistic dynamical models, but current inference frameworks typically assume that reported measurement times are exact. In practice, actual sampling times may deviate from reported times because of sample-handling delays, imper-fect synchronization, or reporting errors. Here, we present a Bayesian framework for parameter inference in ordinary differential equation models that explicitly accounts for uncertainty in measurement times. We formulate latent measurement times as random variables and derive a joint and marginalized posterior. To compute the marginal likelihood efficiently, we augment the original dynamical system with additional state variables that evaluate the required integrals during numerical simulation. This reduces the dimensionality of the estimation problems and allows for efficient and reliable Markov chain Monte Carlo sampling. Across synthetic examples and a published model of carotenoid cleavage in Arabidopsis thaliana, neglecting time uncertainty led to biased estimates and overconfident uncertainty quantification, whereas the proposed marginalized formulation recovered reliable parameter estimates while substantially improving sampling efficiency and scalability. These results identify measurement time uncertainty as an important source of variability in dynamic modeling and establish posterior marginalization as a practical strategy for robust mechanistic inference.

Bacterial communities grow as dynamic populations that respond to their environments. A clinically relevant example is the inactivation of beta-lactam antibiotics by intracellular beta-lactamase in E. faecalis resistant strains. In these populations, resistant bacteria act as antibiotic sinks, detoxifying the environment and allowing sensitive bacteria to survive treatment through a cooperative interaction. In this work, we study strongly coupled planktonic and biofilm populations of mixed sensitive-resistant E. faecalis bacteria under antibiotic stress using fluorescent microscopy. The presence of resistant bacteria in the system benefits both resistant and sensitive cells, leading to mixed planktonic and biofilm populations at super-inhibitory drug concentrations. We show that a beta-lactam antibiotic with or without the addition of a beta-lactam inhibitor can lead to a population inversion effect, characterized by a non-monotonic relation between initial and final fractions of resistant bacteria. The effect is observed in both the planktonic and biofilm populations and is modulated by the total initial cell density. A well-mixed model with competition mediated by resource sharing and cooperation from global degradation of toxins predicts the experimentally observed behavior. These observations suggest underlying population-level mechanisms that are largely independent of biofilm spatial structure.

The basic reproduction number R0 is central to malaria epidemiology, yet it is typically treated as a static quantity derived under memoryless assumptions for mosquito demography. In natural systems, however, mosquito populations are shaped by delayed processes such as larval development and density-dependent feedback, introducing biological memory into vector dynamics. We develop a minimal delay-based framework that incorporates this memory into the Ross-Macdonald model by describing adult mosquito abundance with a retarded differential equation. This formulation induces a time-dependent transmission potential R0(t). Using complex analysis and the argument principle, we derive an explicit stability threshold [Formula], which separates stable from oscillatory transmission regimes. Near this threshold, delayed feedback produces slow relaxation times and sustained transient oscillations, implying that transmission potential may vary intrinsically even in the absence of external forcing. To account for ecological variability, we extend this deterministic condition into a probabilistic framework and define the stability probability as [Formula]. Numerical simulations and global sensitivity analysis show that recruitment and developmental delays are the primary drivers of instability, while adult mortality has a weaker stabilizing effect. These results indicate that malaria interventions may influence not only the magnitude of malaria transmission but also its dynamical stability. By linking delay dynamics, transmission theory, and uncertainty quantification, this framework provides a basis for stability-aware modeling and interpretation of malaria transmission under ecological variability. Author summaryMalaria transmission is often summarized by a single number, R0, treated as a fixed indicator of whether transmission will increase or decline. This assumes mosquito populations respond instantly to environmental conditions. In reality, mosquitoes develop through stages where larval conditions, such as crowding, nutrition, or temperature, affect adult populations only after a delay. This creates biological memory: todays mosquitoes reflect past environments. We show that this memory can fundamentally reshape transmission dynamics. When developmental delays are included, transmission potential is no longer constant but can fluctuate over time, even in stable environments. These fluctuations can persist or amplify depending on the balance between mosquito growth, mortality, and delay. As a result, variability in mosquito abundance or malaria transmission may arise from intrinsic dynamics rather than external drivers alone. Under ecological variability, stability becomes probabilistic, allowing estimation of how likely transmission is to remain stable. Interventions that reduce larval productivity or increase adult mortality may therefore both lower transmission and make it more predictable, improving interpretation and control strategies.

Understanding how biological shape and movement interact with surrounding fluids represents a fundamental challenge at the intersection of biology, physics, and engineering. Fish locomotion exemplifies this challenge: body morphology and swimming kinematics together determine the hydrodynamic forces and flow structures that enable efficient propulsion and maneuverability. Whereas biologists have long sought to connect morphological variation to swimming performance, traditional morphometric approaches provide limited insight into the fluid mechanical consequences of shape differences. Similarly, although computational fluid dynamics can reveal detailed flow physics, simulating hydrodynamics across diverse and dynamic morphologies remains prohibitively expensive for systematic investigation. To bridge this gap, we introduce a data-driven framework that connects fish body shape dynamics to hydro-dynamic performance through compact morphospace parameterization and reduced-order modeling. Using CFD simulations of 15 fish species from the Digital Life Project database (www.digitallife3d.org/3d-model), we generate hydrodynamic datasets capturing the shape-flow relationship. Principal Component Analysis (PCA) extracts four dominant shape parameters from dorsal body profiles, which are then integrated into an Inverse-Design with Dynamic Mode Decomposition (ID-DMD) framework to model the resulting fluid dynamics. The resulting modal analysis suggests that locomotion strategies emerge from specific shape-flow interactions. We further demonstrate the frameworks utility through single- and multi-objective shape optimization, showing how it enables efficient exploration of the morphology-hydrodynamics relationship. This approach offers a novel analysis and design tool for understanding how biological form and motion interact with fluid mechanics, with applications ranging from bio-inspired vehicle development to evolutionary biomechanics.

Temporal gradient sensing is a fundamental capability observed across diverse natural biological systems, contributing to the coordination of their functions. Harnessing this ability is also of significant interest in synthetic biology, particularly for sensing and control applications. In this work, we focus on a biomolecular topology that exemplifies a broader class of signal-differentiating architectures, while introducing a structural variant of it. We examine their behavior under both nominal and non-ideal conditions, accounting for stochastic noise arising from different sources. Our investigation includes scenarios where these topologies operate independently, as well as when embedded within minimal regulatory architectures based on negative as well as positive feedback. We analyze the stability of the resulting macroscopic dynamics--a prerequisite for practical deployment--and quantify stochastic fluctuations in system output, providing comparisons with the corresponding input/unregulated process. Importantly, our results demonstrate that signal differentiation can be effectively implemented in a biomolecular setting without incurring deleterious noise amplification--a major concern in the utilization of derivative action across disciplines.

Gene drives, genetic constructs that can spread deleterious alleles in wild populations, have the potential to address some of the major pressing challenges of the Anthropocene such as invasive species, spread of disease vectors, and agricultural pests. However, responsible and effective deployment of gene drive requires taking into account the complex nature of real-world population connectivity networks. In particular, it is unclear how the topological position of the deployment site affects the spread process and its final outcome. Here we develop a framework for modeling gene drive spread in population connectivity networks, and study the eco-evolutionary dynamics of gene drive spread under complex population structures. We investigated the relationship between the position of the deployment site in the topology of the network and whether the gene drive is eventually lost, fixed, or maintained at an intermediate frequency. We identified network centrality measures of deployment sites that are highly correlated with the outcome of deployment for different gene drive designs and across diverse network topologies. We also show that there is a trade-off between the time-to-fixation and the final outcome, implying that multiple centrality measures of the deployment site would need to be considered when aiming to achieve rapid and successful population control using gene drives.

The ecological and evolutionary dynamics of populations, including viral populations, are known to be jointly shaped by deterministic and stochastic processes. While the impact of stochastic processes has been rigorously explored for viral dynamics at the level of the host population, most dynamic models for acutely-infecting respiratory viral pathogens at the within-host scale remain deterministic in their formulation. While this may be reasonable for identifying key processes shaping their within-host viral population dynamics, recent studies indicate that stochastic processes need to be invoked for understanding patterns of within-host viral evolution. Specifically, several studies have shown that viral allele frequencies can change dramatically over the time course of days in acute infections. Here, we use stochastic dynamic models to explore the role of environmental noise in shaping observed patterns of virus evolution in acute respiratory virus infections. We summarize ways in which environmental stochasticity can be biologically realized in these acute viral infections and describe within-host models that can be implemented to jointly yield viral population dynamics and evolutionary dynamics. We further develop a statistical approach to estimate the extent of environmental noise from observed within-host allele frequency changes. We test this approach on simulated data and apply it to existing influenza A virus and SARS-CoV-2 within-host data. With these applications, we show that environmental stochasticity can parsimoniously reproduce key features of empirically observed allele frequency changes without needing to invoke demographic stochasticity or to adopt Wright-Fisher model formulations with a constant effective population size. Finally, we show that purifying selection and positive selection can both still contribute to within-host viral evolution in the context of a noisy environment, providing theoretical support for studies that have found purifying and positive selection in acutely-infecting respiratory virus populations.

Cytoplasmic mechanical properties are often treated as constant background parameters, yet whether they change systematically during development remains unclear. Here, we directly measured cytoplasmic mechanics during early embryogenesis of Caenorhabditis elegans by establishing active microrheology using micrometer-sized magnetic droplets. Active microrheology revealed a progressive decrease in creep compliance from the 1-cell to the 8-cell stage, indicating a progressive stiffening of the local cytoplasmic environment during development. This decrease persisted even when cytokinesis was inhibited, demonstrating that it cannot be explained solely by geometric changes associated with cell division. Passive microrheology using 40-nm fluorescent beads showed a consistent decrease in probe mobility over development. Together, these results demonstrate that cytoplasmic mechanical properties undergo a gradual, developmentally programmed change during embryogenesis that cannot be explained by cell division-associated geometry alone.

Hydrodynamic synchronization of motile cilia is essential for biological functions such as fluid transport, locomotion, and developmental patterning. It comprises the generation and the response to local flows in complex geometries. Besides their central role in physiology, direct experimental tests of ciliary responses to local flows at cellular length and time scales have remained elusive, largely due to the absence of tools capable of applying controlled, and localized flow stimuli. Here, we introduce programmable, nanometer-thin Ti/Pt microactuators that generate well-defined hydrodynamic forcing at biologically relevant frequencies while operating at biocompatible sub-Volt voltages. This platform is pioneering a controlled local hydrodynamic stimulation of individual motile cilia. We quantify the flow fields and forces produced by single microactuators using particle image velocimetry. Applying local oscillatory flows close to motile cilia of the green alga Chlamydomonas reinhardtii, we probe their dynamic response by quantifying phase-locking between cilia and microactuators. This quantification is aided by combining machine-learning-based image segmentation, oscillator phase reconstruction, and circular statistics. During actuation, we observe signatures of phase-locking: those include a reversible modulation of the fluctuations in phase-difference between cilium and actuator and a systematic shift in ciliary beating frequency. Beyond providing a bio-compatible and precise platform for local hydrodynamic stimulation, our approach establishes an experimental framework for directly testing theories of hydrodynamic synchronization and load adaptation in systems of motile cilia.

Adolescence is an important developmental period during which there are diverse changes in the brain and behavior. Goal-directed behaviors and the component processes underlying those behaviors improve during adolescence, including working memory, response inhibition, and reinforcement learning. At the same time there is substantial pruning of excitatory connections in prefrontal cortex and ongoing myelination of axons. However, psychiatric disorders also become increasingly prevalent in late adolescence and early adulthood. In this study, we develop computational models that suggest a hypothesis for how the ongoing changes in the brain can give rise to the increased prevalence of psychiatric disorders. We show that both myelination and pruning during adolescence lead to attractor landscapes in which strongly encoded memories, driven by three-factor learning rules that modulate Hebbian plasticity, come to dominate the landscape of brain activity, at the expense of weakly encoded memories. Pruning and myelination lead to large, strong attractors which, if they are related to aversive emotions, can drive intrusive thoughts and compulsions in obsessive compulsive disorder, rumination in depression, and aversive memories in post-traumatic stress disorder. The link between pruning, myelination and the emergence of dominant attractors for emotionally salient memories is well supported by the models. The way these effects map onto forebrain circuits requires more work.

Natural populations are often nonlinear and exhibit substantial variability. A central question is how stochasticity interacts with density-dependent regulation to shape population stability. We address this using four long-term time series of Trinidadian guppies and find that their dynamics are well described by a stochastic logistic model with multiplicative environmental noise. The model predicts that stochasticity does not merely add fluctuations around deterministic carrying capacity, but alters the equilibrium structure. Using stochastic bifurcation theory, we show that increasing noise shifts the most-probable population size below the deterministic equilibrium and can push populations closer to a noise-induced bifurcation, even when mean growth rates remain positive. The effects of stochasticity across populations align with known ecological differences among streams, particularly the effects of light level and seasonality. The analysis also identifies populations most sensitive to perturbations, which are not detected by standard early warning indicators. Temporal and spectral analyses further show that intrinsic growth rate governs local recovery, while seasonal variation interacts with density-dependence to shape longer-term population fluctuations. Together, our results show that stochasticity can alter resilience and vulnerability by reshaping ecological stability landscapes.

The classical concept of Critical Community Size (CCS) as formulated by Bartlett defines the minimum host population required for a pathogen to persist endemically without stochastic extinction. While this framework successfully described directly transmitted childhood infections in relatively isolated populations, it is increasingly inadequate for modern urban systems characterized by strong connectivity between cities. Pathogens circulating in highly connected urban networks can repeatedly re-emerge through spatial reintroduction even when local transmission temporarily fades out. In such systems, persistence is inherently probabilistic and influenced simultaneously by population size, environmental suitability, and network connectivity. In this study, we develop a generalization of the CCS concept, the Empirical Persistence Threshold (EPT), and apply it to three of the main arboviruses circulating in Brazil--dengue, chikungunya, and Zika--over the period 2017-2024. The Empirical Persistence Threshold generalizes the classical notion of critical community size by replacing a single deterministic threshold with a probabilistic, datadriven measure. Instead of asking for the minimum population at which persistence is guaranteed, EPT characterizes the lower tail of the population distribution among municipalities that empirically sustain transmission. Using weekly incidence data from thousands of municipalities, we transform temporal incidence series into binary sequences describing the presence or absence of reported transmission. For each municipality, we characterize persistence through the empirical distribution of run lengths of consecutive weeks with reported cases. Distances between run-length distributions are computed using the Wasserstein-1 metric, allowing a geometrically meaningful comparison between persistence profiles, and municipalities are grouped into epidemiological regimes using hierarchical clustering methods. Across all three arboviruses, we identify two robust regimes: one exhibiting sporadic and recurrent epidemic transmission, and the other exhibiting sustained persistent transmission. We then estimate the population scales associated with each persistence regime. The analysis is further extended to evaluate how persistence thresholds vary across climate regimes (Koppen classification) and urban hierarchy levels (REGIC). This framework allows the estimation of probabilistic persistence thresholds analogous to CCS, but adapted to connected urban systems. We define the Empirical Persistence Threshold as lower quantiles of the population distribution among municipalities in the persistent regime, and additionally estimate persistence thresholds based on regime membership probabilities. Results reveal strong interactions between population size, climate, and urban connectivity. Dengue exhibits the lowest persistence thresholds, Zika intermediate thresholds, and chikungunya the highest thresholds. These findings demonstrate that pathogen persistence in modern urban systems cannot be described by a single deterministic population threshold. Instead, persistence emerges from the joint effects of demographic scale, environmental suitability, and network position within metapopulation systems. Author SummaryInfectious diseases often require a minimum population size to persist locally, a concept known as the critical community size (CCS). This idea was developed for relatively isolated populations, but modern cities form highly connected networks where diseases can repeatedly reappear even after local transmission disappears. In this study, we introduce the Empirical Persistence Threshold (EPT), a data-driven approach that replaces the idea of a single fixed threshold with a probabilistic description of persistence. Instead of focusing on case counts, we analyze how long transmission persists over time in each municipality. Using weekly data for dengue, chikungunya, and Zika across Brazil from 2017 to 2024, we identify distinct patterns of transmission persistence and estimate the population levels associated with sustained transmission. We also examine how these thresholds vary with climate and urban structure. Our results show that persistence depends not only on population size, but also on environmental conditions and the position of cities within the urban network.

When only some hosts are protected from disease vectors, disease spread may be inhibited through a net reduction in vector visits or amplified as vectors redirect their attention to unprotected hosts. Two factors that determine which outcome prevails are host microbiota that alter vector host-seeking behavior and natural enemies that redistribute or suppress vector populations. Because both shape the frequency and distribution of vector visits, they are essential for understanding how individual-level protection scales to population-level disease dynamics. Yet, how these processes interact across scales remains poorly understood. Pea aphids are major virus vectors in pea crops and are commonly managed using parasitoid wasps. Recent evidence suggests that epiphytic bacteria in the genus Pseudomonas can also repel or kill pea aphids, yet whether Pseudomonas complements or undermines parasitoid-based vector control remains unknown. We used a mathematical model to show when and why Pseudomonas complements versus undermines biocontrol of aphid-vectored virus outbreaks. The effect of Pseudomonas on virus outbreaks depends most strongly on how successful parasitoids are at tracking aphid densities: When parasitoids effectively track aphids, Pseudomonas inhibits virus outbreaks by reducing aphid densities. With poor parasitoid tracking of aphids, Pseudomonas-induced aphid mortality generates spatial variability in aphid densities that slows parasitoid population growth. The net result is amplified crowding in plants not protected by Pseudomonas, increasing winged aphid production and accelerating viral spread. Counterintuitively, the more effective Pseudomonas is at killing aphids, the more strongly it generates spatial variability and promotes virus spread. The only other factor that can change the direction of Pseudomonas effects on virus outbreaks is whether the virus starts on a Pseudomonas-protected plant, which can cause Pseudomonas to inhibit virus outbreaks when it would otherwise promote them. Our results show how community and spatial context dictate whether microbiota protective to individual hosts will accelerate viral outbreaks.

Apohas Liquid State Intelligence Platform (LSIP) records ellipsometric waveforms from injections depositing sub-microgram quantities of antibody drop-by-drop onto a liquid reservoir. We previously showed that a behavioural feature extracted from the waveforms, VIBE1, identified antibodies carrying multiple biophysical liabilities in an industrial dataset of 71 monoclonal antibodies, and enriched for clinical failure across a larger dataset of 235 therapeutic antibodies [1]. Here, we use an auxiliary coalescence-sensor channel to decode VIBE1 by separating the coalescence event from its propagation through the substrate. The pertitration drop-to-drop standard deviation of pinch-off time,{sigma}{tau} , explains most of VIBE1s variance across the dataset (R2 = 0.92, n = 1182). High-speed imaging at 10,000 frames per second reveals that all imaged drops initially thin at the same Newtonian capillary-inertial rate while the neck remains wide. In drops from certain antibodies, the thinning bridge then decelerates as internal strain builds in the narrowing neck. This elasto-capillary stiffening response has a timescale{lambda} that decreases as pinch-off time{tau} i increases across the imaged set.{sigma}{tau} is therefore a readout of the antibodys propensity to undergo a transient gel-like stiffening response during coalescence, and that variability is what VIBE1 captures. The signal is concentration dependent, and absent in bovine serum albumin (BSA) tested at up to an order of magnitude higher molarity than the antibodies, despite BSA being a strongly surface-active globular protein. The instrument is configured so that complex behaviours of this kind appear in its recorded waveforms; the gel-like coalescence response we identify here is one such phenomenon.

Cooperative and competitive interactions among individuals harvesting resources can shape environmental states, such as prey abundance. In turn, environmental conditions feed back to influence strategic interactions. Eco-evolutionary game theory studies how these feedbacks shape the co-evolution of behavior and environment. Existing models typically assume deterministic, noise-free environmental dynamics. However, real environments are inherently stochastic, for example due to finite resources, and noise can qualitatively alter social outcomes. Here, we incorporate stochastic environmental dynamics into eco-evolutionary game theory. When environmental change is slow relative to strategy updates, we show that behavior reflects a mixture of the games associated with low and high environmental states, often yielding outcomes qualitatively distinct from deterministic predictions. In particular, environmental stochasticity can eliminate bistability and enforce dominance of a single behavior. When environmental dynamics are faster, populations have less opportunity to track fluctuations, and behavior converges toward strategies that are optimal on average. Stochasticity can even causes persistent oscillations in the tragedy of commons, in regimes where classical models predict stability. Our framework provides a tractable approach for analyzing social behavior linked to environmental dynamics how noise shapes long-term eco-evolutionary outcomes.

Take-off is a fast and energy-efficient strategy for bipedal animals, such as birds, to achieve rapid movement; however, how muscle physiology scales to govern this universal behavior remains unresolved. Research in other species physiologies is not readily applicable. As a result, important questions, whether theropod dinosaurs such as Tyrannosaurus rex were capable of jumping, remain unanswered. In this article, we coupled Lagrangian dynamics with Hills muscle equations and developed new experimental methods to quantify joint rotational stiffness and damping, thereby enabling a systematic description of lower-limb mechanics. The approach establishes a novel kinetic framework that links muscle contractile properties to lower-limb performance without invoking control optimization. Animal observations and tabletop mechanisms validate the framework. The mechanics model reveals that the take-off time of about 0.1 s across body masses of 0.003 to 90 kg is achievable, as heavier birds generate proportionally higher reaction forces. Additionally, Tyrannosaurus rex should be capable of jumping, based on the available physiology data. Beyond evolutionary insights, our framework provides a new methodology for analyzing the mechanical properties of biological joints and informing the design of scalable bio-inspired robots.

Metabolic networks are typically viewed as homeostatic systems that stabilize flux, energy charge, redox balance, and metabolite availability under perturbation. However, it remains unclear whether the same feedback architectures that support metabolic robustness can also generate learning-like, experience-dependent adaptation. Here, we develop a coarse-grained dynamical model of mammalian energy metabolism to test whether prior perturbation can improve future metabolic responses. The model represents core glucose, glutamine, fatty acid, and oxidative phosphorylation pathways as coupled ordinary differential equations with Michaelis-Menten-type fluxes, product-inhibition feedback, adaptive enzyme-capacity regulation, and explicit ATP costs for enzyme adjustment. Rather than aiming to reproduce quantitative fluxes for a specific cell type, the framework is designed to expose how metabolic feedback, regulatory cost, repeated perturbation, and environmental variability interact. We use this model to ask whether adaptive enzyme regulation enables improved recovery after repeated challenges, whether such effects depend on energetic control costs, and whether environmental variability broadens or constrains the set of reachable adaptive states. This approach provides a tractable way to investigate how homeostatic metabolic regulation may give rise to experience-dependent metabolic plasticity.

Inferring gene regulation from time-course expression profiles is essential for understanding how cells transition between states during development, differentiation, and disease progression. Existing approaches often model expression dynamics with ordinary differential equations (ODEs). However, due to the computational complexity of directly solving these ODE models, most methods rely on finite-difference approximations of temporal derivatives, which can amplify measurement noise, introduce discretization bias, and lead to unstable or biased parameter estimates. To fill this gap, we develop the first computational method to directly learn a linear ODE model for gene regulation inference without relying on finite-difference approximations. We first formulate an optimization problem that directly exploits the closed-form solution of the linear ODE system. We then solve this problem via gradient descent, deriving analytical gradients with respect to the model parameters; these gradients involve matrix exponentials and integrals, which are challenging to directly compute. To make the computation efficient, we further use high-order Taylor approximations of the gradients whose truncation error is on the order of machine precision. In addition, we establish theoretical results demonstrating an inherent, non-vanishing gap between our exact solution and solutions derived from finite-difference approximations, which underscores the theoretical advantages of our approach. Finally, we demonstrate that our method consistently outperforms competing approaches on both simulated data and real-world scRNA-seq datasets in terms of AUROC. Our source codes can be accessed here: https://github.com/EJIUB/ExactLinearODE

1Although resources are typically distributed continuously in space, species distributions often organize into discrete clusters. In his seminal paper [36], Turing demonstrated that such clusters can spontaneously arise in population densities, even when populations evolve in environments with continuously varying conditions. This phenomenon is known as Turing instability. In this work, we focus on two models grounded in population dynamics: a one-dimensional model based on the nonlocal Fisher-KPP equation, and a two-dimensional model involving an environmental gradient. We show that phenotypic clusters (sometimes referred to as "species") emerge in these models. We prove that they do not emerge because of Turing instability, but because of stochasticity, and that they disappear when stochasticity is reduced. First, for both models, we start our simulations with initial populations uniformly distributed in the state space. We show that phenotypic clusters quickly emerge and that the distances between them depend on the population size, that is, on the degree of stochasticity. Next, we start from already clearly defined phenotypic clusters. We identify three regimes in the connection between population size, the initial distances between clusters, and the distances between clusters at equilibrium. Last, on the two-dimensional model, we relax the hypothesis of complete clonality by varying the effective recombination rate, explore its effect on phenotypic clustering, and show that phenotypic clustering decays drastically with slight recombination.